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Techniques for determination of very precise orbits for satellites of the Global Posi- 
tioning System (GPS) are currently being studied and demonstrated . These techniques 
can be used to make cm-accurate measurements of station locations relative to the geo- 
center, monitor earth orientation over timescales of hours , and provide tropospheric and 
clock delay calibrations during observations made with deep space radio antennas at sites 
where the GPS receivers have been collocated. For high-earth orbiters, meter-level knowl- 
edge of position will be available from GPS, while at low altitudes, sub-decimeter accu- 
racy will be possible. Estimation of satellite orbits and other parameters such as ground 
station positions is carried out with a multi-satellite batch sequential pseudo-epoch state 
process noise filter. Both square-root information filtering (SRIF) and UD- factorized 
covariance filtering formulations are implemented in the software. A Bierman-modified 
Rauch-Tung-Striebal (BRTS) smoother runs in conjunction with the SRIF and UD fil- 
ters to compute smoothed estimates and covariances. The filtering algorithms have been 
arranged to take advantage of sparse matrices and other characteristics of the GPS mea- 
surement scenarios. The filter includes unique error evaluation capabilities to assess 
effects from mismodeling. Process noise plays a key role in the orbit determination for 
stochastic behavior of transmitter/receiver clocks, atmosphere-induced delay fluctua- 
tions, and unmodeled satellite accelerations. The efficiency and accuracy of the SRIF 
and UD filter formulations are compared for GPS processing under a variety of condi- 
tions. With data from recent GPS experiments using the seven satellites currently in orbit , 
continental ground baselines have been measured with GPS and with VLBI which agree to 
within 2.5 cm over distances of 2000 km, corresponding to a relative baseline accuracy of 
better than 1.5 parts in 10 8 . 


I. Introduction 

The Global Positioning System (GPS) will consist of at least 
21 satellites launched by the United States Department of 
Defense equally spaced in six orbit planes at about 20,000 km 
altitude. The satellite constellation, currently consisting of 


seven operating satellites, is expected to be complete by the 
early 1990s. Many types of scientific applications are uniquely 
suited to the precise positioning capabilities provided by GPS. 
The relatively high precision, low cost, mobility, and conven- 
ience of GPS receivers make GPS-based positioning attractive. 
Dense GPS ground networks are already operating in North 
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America, South America, Europe, and Japan to monitor cm- 
level crustal motions in geologically active regions. Initial results 
from GPS experiments since 1985 are encouraging and suggest 
that accuracies may be obtained equaling those available from 
other generally more restrictive techniques, such as Very Long 
Baseline Interferometry (VLBI) or satellite laser ranging 
(SLR). When the GPS constellation is complete, a worldwide 
ground tracking network equipped with advanced GPS receivers 
will enable sub-decimeter low-earth orbiter position accuracies 
to be achieved for satellites such as TOPEX/Poseidon fl] and 
the Earth Observing System (Eos) platforms [2, 3) . This world- 
wide tracking network will also be used to monitor earth 
orientation over timescales of less than one day [4] , improve 
relative station positions at the cm level, determine absolute 
station positions relative to the geocenter with accuracy of 
better than 5 cm [5] , enable global time transfer at the nano- 
second level, and provide precise calibrations of path delays 
from the ionosphere and troposphere. The improved knowl- 
edge of earth orientation, station positions, time transfer, and 
media calibration will be used at NASA’s Deep Space Network 
(DSN) tracking stations in support of navigation for planetary 
exploration missions. The GPS-based tracking network will 
also provide a meter-level orbit determination capability at the 
DSN for certain high -earth orbiters [6] . 

In this article, an overview of GPS signal structure and posi- 
tioning techniques will be presented. A section describing the 
filter/smoother software developed and used at the Jet Propul- 
sion Laboratory (JPL) for GPS orbit determination and param- 
eter estimation will highlight some of the filtering techniques 
and algorithms which are especially useful in GPS-based least- 
squares parameter estimation. Recent high-accuracy GPS orbit 
and ground baseline results will also be presented to demon- 
strate the feasibility and potential of high-accuracy GPS-based 
navigation and positioning applications. 


II. GPS Signal Structure and Positioning 
Techniques 

The Global Positioning System is designed so that typically 
four to eight GPS satellites are visible simultaneously at any 
time from most locations in the world. The GPS satellites 
transmit down towards the Earth carrier signals at two L-band 
frequencies (1.227 and 1.575 GHz), which are modulated by 
a pseudorandom noise code, the P-code (Fig. 1). The two fre- 
quencies enable the user to remove most of the signal delay 
originating in the ionosphere. When four satellites are in view, 
the user equipped with a receiver which can receive the GPS 
P-code has enough information to solve for the user position 
and the user clock offset from GPS time. This is the simplest 
and most basic GPS positioning technique, often called direct 
GPS tracking (Fig. 2a). A second code, the coarse acquisition 


(C/A) code meant for acquiring and locking onto the GPS 
signal, can also be used for user positioning, but with some- 
what degraded accuracy compared to the P-code. Part of the 
navigation message broadcast down on the P-code includes the 
ephemeris for each satellite and the clock offset of each satel- 
lite from the reference GPS time. The accuracy of the broad- 
cast ephemeris (about 10-20 m) and the clock information 
currently limit user positioning accuracy to about 10-15 m 
with real-time techniques using the P-code. The P-code observ- 
able is often called a pseudorange since it measures the range 
to the satellite but includes offsets between the transmitter 
and receiver clocks. 

To improve positioning accuracy through cancellation of 
clock errors and partial cancellation of ephemeris errors, 
differential GPS tracking techniques are employed (Fig. 2b). 
Instantaneous P-code positioning accuracy of 1-3 m is possible 
with differential GPS tracking. For certain military applica- 
tions as well as for some real-time applications such as an earth 
orbiter docking with the space station, accuracies obtainable 
with instantaneous direct or differential GPS tracking are 
sufficient. 

There are several strategies for reducing the GPS error bud- 
get further for high-accuracy applications, which include cm- 
level station determination, meter-level high-earth orbit estima- 
tion, cm-level earth orientation monitoring, and sub-nanosec 
time transfer. The high-precision strategies include use of 
precise GPS carrier phase tracking and averaging over several 
hours or more, use of dynamical information to reduce the 
number of degrees of freedom, and estimation of improved 
GPS orbits and clocks along with determination of other esti- 
mated parameters from the tracking data. The GPS carrier 
phase provides a very precise measure of range change , but is 
ambiguous in absolute range to an integer number of carrier 
wavelengths. Measurement noise with modern GPS receivers 
over five -minute intervals is typically several mm with the dual 
frequency combined carrier signal (to remove ionospheric 
effects). If the satellites are tracked for several hours, the sig- 
natures in the data enable estimation of the integer cycle ambi- 
guities along with the other parameters. A more powerful 
approach is to process the P-code pseudorange, which can be 
one or two (or more) orders of magnitude noisier but provides 
a measure of absolute range, along with the carrier phase. This 
constrains the clocks and carrier phase range ambiguities so 
that the solutions are strengthened considerably. 

Most commercial receivers provide pseudorange with mea- 
surement scatter of 60-200 cm over five -minute averaging 
intervals. Most of this scatter is due to multipath character- 
istics of the ground site and antenna. The Rogue receiver was 
designed at JPL with the goal of substantially reducing pseu- 
dorange measurement scatter. The Rogue, which features 
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digital baseband electronics, can track up to eight satellites 
simultaneously with dual-frequency carrier phase and pseudo- 
range. A new antenna design shows considerable promise in 
reducing P-code multipath by an order of magnitude (Fig. 3). 
This design includes a drooped cross-dipole antenna with a 
non-absorbing backplane which has concentric choke rings 
machined from an aluminum disk spaced to maximally elimi- 
nate multipath signals. When GPS equipment such as the 
Rogue receiver and its new antenna design are routinely used 
in the field, a substantial improvement in tracking perfor- 
mance and system accuracy is expected. 

III. Formulation of the Least-Squares 
Problem 

As discussed above, to improve GPS-based tracking perfor- 
mance. GPS carrier phase and pseudorange can be received 
over a period of hours and parameters then estimated using 
dynamical model information. At JPL. the GPS orbit determi- 
nation problem is reduced to a system of linear equations 
which can be solved with a batch-sequential process noise fil- 
ter. This approach requires a moderately accurate nominal 
model for the GPS orbits and station locations so that only the 
linear terms in an expansion for deviations from this nominal 
model need be retained. It has been found that initial GPS 
epoch states should be accurate to about 1 km and receiver 
locations should be known to within a few hundred meters. 1 
The initial GPS epoch states are numerically integrated with a 
multi-step Adams method 2 in the J2000 inertial reference 
frame [7] . Variational partial derivatives are computed relat- 
ing the change in satellite position and velocity with respect to 
changes in the initial epoch states and with respect to force 
parameters. The Earth's gravity field is expressed in terms of a 
spherical harmonic expansion, and the gravitational effects o i 
the sun, moon, and planets are represented as due to point 
masses. The JPL software includes the ROCK4 [8-10] GPS 
solar radiation pressure model, allows for estimation of arbi- 
trary unmodeled accelerations on the satellite, and includes 
impulsive motor burns when needed to model GPS maneuvers. 
Typically the broadcast ephemeris or a more precise ephemeris 
produced by the Naval Surface Weapons Center (NSWC) is 
used for the nominal trajectory. 

Precise earth models are used both to model the measure- 
ments and to compute the partials for measurements with 


*If the nominal model is not accurate enough, the output from the fil- 
ter can be used as the new nominal model and the solution can be 
iterated. 

2 F. T. Krough, "‘Changing Stepsize in the Integration of Differential 

Equations Using Modified Differences,” JPL Section 914 TM No. 312 

(internal document), March 20, 1973. 


respect to model parameters [11] as a function of time. The 
algorithms used are based on ones developed for Very Long 
Baseline Interferometry (VLB1). an extremely accurate astro- 
metric technique which uses radio astronomical measurements 
of fixed point sources (quasars) [12] for deep space naviga- 
tion and determination of ground baselines and earth orienta- 
tion. The models include UT1-UTC. polar motion, nutation, 
precession, solid earth tides, ocean tidal loading, general rela- 
tivistic clock corrections, and the Lanyi tropospheric delay 
mapping function [13], which relates the tropospheric delay 
at various elevations to its value at the zenith. Through the 
chain rule, the partials are referred to the GPS epoch states. 
The GIPSY (GPS Inferred Positioning SYstem) software is 
used for GPS data processing. The OASIS (Orbit Analysis and 
Simulation Software) software, developed at JPL 3 ' 4 for co- 
variance analysis and simulations, shares most programs with 
GIPSY but has a streamlined, simplified models module. 

The difference, z, between each observation and its com- 
puted predicted value from the nominal model is calculated. 
The measurement equation is 

z = Ax + Ap + Ay + p (1) 

x p r y 

where A is a row from the matrix of measurement partials, x is 
the satellite state vector (usually three position and three veloc- 
ity components), p is the vector of process noise parameters, 
y is a vector of bias (constant) parameters, and v is a zero 
mean white noise with iv 2 ) equal to the measurement noise 
variance. The OASIS/GIPSY filter is a multi-satellite program, 
so each satellite in turn has its associated state vector = 
(x T .p T ,y r ) T , and these are ordered one after the other in the 
total state vector X, with common parameters (such as station 
locations) placed at the end. The measurements are combined 
and an optimal least-squares solution for the estimated param- 
eters determined in the filter. The filter divides the measure- 
ments into finite time intervals known as batches. Within each 
batch, all process noise parameters-parameters which are 
modeled as stochastically time-varying-are assumed to be 
piecewise constant. After the filter processes the measure- 
ments in a given batch, to update the estimates and covariance 
for the parameters a time update is performed to fold in the 
effects of process noise. Measurements are then processed in 
the next batch, and so on. After filtering is completed, a 
smoother works recursively backwards in time to update op- 


3 S. C. Wu, W. I. Bertiger, J. S. Border, S. M. Lichten, R. E. Sunseri, 
B. G. Williams, and J. T. Wu, OASIS User's Guide, V 1.0, JPL D-3138 
(internal document), April 1, 1986. 

4 S. C. Wu, W. I. Bertiger, J. S. Border, S. M Lichten, R. F Sunseri, 
B. G. Williams, P. J. Wolff, and J. T. Wu, OASIS Mathematical Descrip- 
tion, V. 1.0, JPL D-3 139 (internal document), April 1, 1986. 
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timally computed estimates and covariances when the pres- 
ence of process noise requires it. 

Both the OASIS and GIPSY software use the same filter 
for estimate/covariance computation. The filter offers the user 
the choice of two different mechanizations: the square-root 
information filter (SRIF), or the UD-factorized filter. Both 
mechanizations are based on the algorithms developed in [14] 
by Bierman. The smoother is formulated in terms of UD back- 
wards pass recursions and works with output from either the 
SRIF or UD filter. Both the filter and the smoother make 
heavy use of the Bierman Estimation Subroutine Library 
(ESL) [15] . 

A. UD Filter Formulation for Measurement 
Processing 

The conventional Kalman filter has a measurement update 
mechanization based on the covariance matrix , P. Let the esti- 
mate vector and covariance matrix be X and P, with the con- 
vention that ~~ denotes quantities before measurement updating 
and ~ denotes quantities after the measurement update. Then 

X = X + KO-AX) 

= X + (G/a) (z - AX) (2) 

P = P-(l/a)GG T 

where A is a row vector of measurement partials, G is the 
normalized Kalman gain vector, which is related to the un- 
normalized Kalman filter gain, K, by 

G = PA t = aK (3) 

and a = apa t + r is the innovations, or prediction residual 
variance. The measurement variance is r. It is assumed here 
that the vector of observation errors, v {y is described by 
E{Vi) = 0,E(vv T ) =Pj,,for / = 1 , . . . , ra,andP^ = Diag(rl , . . . , 
r m ) for the measurement variances. In the case of correlated 
measurement noise, whitening procedures can be used to 
uncouple the measurements and yield diagonal P v [14] , so the 
formulas presented here are valid without loss of generality. 

The conventional Kalman measurement update (Eq. 2) has 
been shown to be sensitive to computer roundoff with some- 
times catastrophic loss of accuracy from imperfect cancella- 
tion when positive quantities are differenced. Conventional 
Kalman algorithms have also exhibited sensitivity to numerical 
ill-conditioning [14, 16] . A comprehensive comparison of 
conventional, stabilized, and factorized formulations of the 
Kalman filter for a spacecraft navigation problem [16] demon- 
strated convincingly the superior numerical characteristics of 


the factorized Kalman filter approach. Ill-conditioning with 
unfactorized Kalman filters is sometimes attributed to large 
a priori state uncertainties and relatively small data covariances. 
These very attributes, however, typify many GPS earth orbiter 
problems relying on the GPS carrier phase data type: the car- 
rier phase data variance is on the order of 10“ 5 m 2 while the 
a priori position covariance could be between 1 and 10 6 m 2 . 
Another filter characteristic associated with ill-conditioning is 
a very low level of process noise. However, low process noise 
levels are common in high-precision GPS applications when 
stochastic models are used to model slowly varying tropo- 
spheric delay variations, GPS force parameters, and stable 
hydrogen maser station clocks. Near-singular covariance 
matrices are sometimes encountered in GPS orbit filtering 
because of different receiver characteristics. At some sites, 
GPS receivers were operated side by side in experiments which 
took place in 1985 where one receiver was capable of pro- 
ducing GPS carrier phase and P-code pseudorange and the 
other could only produce carrier phase data (codeless). The 
clock offset for the code receiver could be determined much 
more precisely than for the phase-only receiver: the clock 
variances were typically 2.5 X 10“ 7 /isec 2 for the code re- 
ceivers while codeless receiver clock variances are typically 
10 4 jusec 2 . Furthermore, since these clock parameters are 
often modeled with process noise, their variances vary con- 
siderably in the course of a run as satellites move through 
different geometries with good and bad observability. With 
GPS scenarios, severe numerical difficulties were sometimes 
encountered when mapping orbits and baselines to various 
coordinate systems when using the unfactored covariance 
matrix; with factored or square-root form, however, these 
mappings were stabilized and no numerical problems have 
been noted. Experience at JPL with deep space tracking esti- 
mation software has also resulted in complete conversion to 
factorized (or square-root) Kalman filter mechanizations in 
the orbit determination program used for missions involving 
such spacecraft as the Voyager, Magellan, and Galileo inter- 
planetary probes. With efficient algorithms [14-16], there 
is only an insignificant penalty in computation speed when 
using factorized instead of conventional Kalman filters. 

An extremely stable factorized version of the Kalman mea- 
surement update equations exists for the U and D factors of 
the covariance matrix, where U is an upper triangular matrix 
and D is a diagonal matrix: 

P = UDU t (4) 

The method used is an upper triangular UD Cholesky 
square-root-free factorization with the geometric form reduc- 
tion described in [14]. In order to save storage, U is vector 
stored, and the factorization is chosen so that U has unity on 
its diagonals and D can be stored in those locations. It is con- 
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venient to store the estimates, X, in the last column appended 
on the end of the UD matrix. 

With the UD factorization, the measurement update can be 
expressed as 

P = UDU t = U(D - cW T )U T (5) 

where 

V = DF 
F = (AU) t 

The scalar c - 1/a = 1/(7" + V^F) is computed from the 
prior covariance and the current measurement noise. Thus, the 
normalized Kalman gain can be expressed in terms of the UD 
factors as 

G = UDU t A t = UV (6) 

The measurement update algorithm used in the OASIS/ 
GIPSY filter is the Bierman-Carlson-Gentleman update [14, 
16, 17]. This algorithm avoids explicit differencing when 
computing updated diagonals because explicit differencing for 
these terms makes the update prone to loss of accuracy due to 
imperfect cancellation, leading to negative diagonals in extreme 
cases. 

The UD filter processes measurements one at a time. Betore 
updating, each measurement is tested for consistency with the 
prior estimates and covariance. If the measurement fails the 
test, it is considered to be an outlier and not processed. The 
prediction residual must satisfy the criterion 

(z-AX) 2 < TSTot (7) 

where TST is the level of acceptance and a is the innovations 
variance defined above. For example. TST = 9 is a 3 -o test. It 
Eq. 7 is satisfied, the estimates and covariance are updated. 

B. SRIF Formulation for Measurement 
Processing 

The SRIF data equations are 

z = RX + v (8) 

where 

P = R-'r t 


and estimates are retrieved as 

X = R -1 z 

The square-root information array [R:z] is analogous to 
the UD-factored covariance with estimates in the last column 
of the UD matrix. The square-root information matrix, R, like 
the U matrix, is upper triangular and vector stored. Note that 
for the SRIF, the noise v is zero mean as before but has unit 
covariance, so z and v have been normalized by the measure- 
ment noise. 5 The SRIF data processing algorithm in the 
OASIS/GIPSY software follows [14] and [18] in applying 
orthogonal Householder transformations, T H , to an aug- 
mented information array 


1 

/ 

dfL 


R z 

M 


0 e 


The measurement partials (normalized) are in A, and e con- 
tains information about the post-fit residuals. The matrix A 
ordinarily holds a subset of the data, since on some computers, 
available memory or disk space limits how many measure- 
ments can be processed at once. When process noise is in- 
cluded, the measurements to be processed together in a buffer 
must lie in the same time batch. The ith row in A contains the 
/th measurement in this measurement buffer, so when the / th 
elementary Householder transformation is applied, only the 
/th row of R and all columns to the right of column i in A are 
affected. The Householder transformation, as implemented in 
the filter, maintains the upper triangularity of the R matrix 
[14, 18]. 

The effect of the size of the measurement buffer on the 
speed of the SRIF on a DEC Microvax II is shown in Fig. 4. 
Beyond a buffer of about 50 measurements, little increase in 
speed is seen. However, for very small buffer sizes, the House- 
holder transformation becomes inefficient. The reason for this 
is that, based on the operation counts in [14] , for additions 
and multiplications the ratio of (overhead)/ (total processing 
time) in processing each measurement buffer with M b measure- 
ments is 3/(3+2 M b ) and \j(\+2M b ), respectively. 

There is also an additional amount of overhead from the N 
square roots, N divisions, 2N additions, and N multiplications 
needed for each buffer being processed, where N is the number 
of parameters estimated. The penalty from these extra opera- 
tions depends somewhat on the relative CPU speed of these 


5 To distinguish these normalized quantities Irom the unnormalized 
ones from the discussion of the UD filter, we represent them without 
italics. 
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operations. The case shown in Fig. 4 corresponds to 3666 mea- 
surements with 192 estimated parameters for a SRIF running 
on a DEC Microvax, where square roots and divisions are 
about 8 and 1.5 times as expensive as additions or multiplica- 
tions, respectively. The composite ratio of overhead to process- 
ing time was calculated and would be 60 percent for this run 
with M b = 1, and less than 1 percent for M b = 100, in rough 
agreement with the curve in Fig. 4. 

C. Process Noise 

Stochastic processes in the GPS orbit software are assumed 
to be piecewise constant over a specified batch interval. Cur- 
rently, first-order Gauss-Markov random processes models are 
used. At the end of a batch, a process noise time update adds 
noise to the covariance matrix and thus causes the time-varying 
behavior of the stochastic parameters. The process noise time 
update for the /th batch maps the estimates and covariance for 
the stochastic parameters into the batch y+1 : 


x . = x. + <t> (/ )p 

y+l } p' J 

Py+i = M /P/ + w / (10) 

p /+1 = <M\<t» T + Q 

where 
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and, as in Eq. (1), the estimated parameters are partitioned 
into satellite states, process noise parameters, and bias (con- 
stant ) parameters 


x 


X = 


P 

y 


do 


M is a diagonal process noise mapping matrix. The process 
noise w ; is a random process with zero mean and for the /th 
process noise parameter. 


E{w>..w. v ) - q.h , (12) 

v ij ik y ^ ij jk ' ' 

with the covariance matrix Q zero except for diagonal ele- 
ments q i in the j th batch, and b jk is the Kroneeker delta 
function. The diagonal entries of M are given by 


m a = ex P Hf, + , -t.)/T..} (13) 


where is the start time for the /th batch and r f . is the time 
constant for the /th stochastic parameter at the /th batch. The 
corresponding diagonal entry in the matrix Q is 

4 , = (1 -m 2 )o 2 ( 14 ) 

where u iss , the steady -state sigma for the /th stochastic param- 
eter, is the noise level which would be reached if the system 
were left undisturbed for a time much greater than r. 6 For 
white process noise, r = 0. m = 0, and as can be seen in Eq. ( 10), 
the filter resets the covariance for the process noise parameters 
at the end of each batch, including zeroing of off-diagonal 
terms and putting in q for the variance on the diagonal. The 
opposite limiting case is the random walk: here both a and r 
are unbounded (r — °°) and a steady-state is never reached. For 
the random walk, M is equal to the identity matrix, and it is 
the rate of change of the process noise covariance, q - dqidt = 
Ar//A/, which characterizes the process, where At is the batch 
size and Aq is the amount of noise added per batch. The Allan 
variance [19]. o\ (Ar), which is often used to characterize 
clock and atmospheric fluctuations [12], can be defined in 
terms of q: 


o 2 a (Ar) - q/At (random walk) (15) 


Thus a random walk process has a slope of -1 for the log-log 
relation between the Allan variance and At. 

The matrix 4> (7). represents the deterministic portion of 
the time update and is nonzero only when there is a dynamic 
coupling between process noise parameters and satellite states. 
M contains the stochastic portion of the time update, x repre- 
sents the pseudo-epoch state of the spacecraft [14] . The cur- 
rent state , x(r ). would be mapped according to 
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The pseudo-epoch state . x ; . is defined relative to the 
current state, x( t ), from 


^Although <jf X and r can vary with time, the subscript j has been left 
off l or simplicity of notation. 
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x(t.) = VVVV Vvv )y 
*/ = *- x 1 (r r ‘ 0 )[x(t J )-* y (‘ r 'o)y] 


(17) 


With this definition, if the pseudo-epoch state is used in the 
filter instead of the current state, the time update of Eq. (10) 
results: 
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(18) 


The mapping matrix (/) used in the pseudo-epoch state 
filter is calculated from the variational partials: 


*//)■ [■Mv.-V 1 ' 1 VW'o> 


3x 




7+i 




P/ = o 


(19) 


D. Smoothing 

The OASIS/GIPSY smoother incorporates Bierman’s modi- 
fication of the Rauch-Tung-Streibel (RTS) smoother [21 ], 
known as the BRTS, including the new reformulation of 
Bierman’s original UD smoother [22, 23] . The smoother 
algorithm makes use of a decomposition of the linear model 
dynamical equations so that rank-1 and rank-2 matrix modifi- 
cations can be substituted for more complex (and costly) 
covariance updates. The smoothing algorithm uses the UD 
formulation, is based on numerically stable Givens reflections, 
and runs “backwards” in time using a recursive method for 
generating smoothed estimates and smoothed UD factors. The 
last (terminal) filter UD or SRIF matrix is used to initialize the 
smoother, so in the case of the SRIF, the terminal square-root 
information array is inverted to UD form. The smoother has 
considerable flexibility in the handling of singular and near- 
singular covariance matrices. As with the deterministic and 
stochastic filter time updates, the OASIS/GIPSY implementa- 
tion of Bierman’s smoothing algorithm takes advantage of the 
sparse matrix structure which typifies GPS scenarios and 
orders parameters to minimize CPU cycles, memory require- 
ments, and disk storage. 

E. Evaluation of Filter Mismodeling 


The use of the pseudo-epoch state formulation saves con- 
siderable computation in the filter and smoother since the UD 
time update is accomplished with sparse matrix multiply rou- 
tines. When no dynamic stochastic parameters are present, the 
pseudo-epoch state is the same as the true epoch state. Because 
of the pseudo-epoch state formulation, most of the mapping 
matrix is filled with zeros or ones (see Eq. 18). The state vec- 
tor X is arranged so that the satellite state for each satellite 
and the associated stochastic parameters are grouped at the 
top. A subroutine loop performs the sparse matrix multiply to 
update the U matrix and estimates for each spacecraft in turn, 
making efficient use of knowledge of the implicit zeros and 
ones for the deterministic mapping. After the deterministic 
4> p mappings for all the satellite are finished, the stochastic 
M mappings are performed to update the process noise states. 
The UD stochastic update is accomplished with the Bierman- 
Thornton one-at-a-time update [14, 17, 20] . The deterministic 
portion of the SRIF time update also takes advantage of the 
sparse mapping matrix structure. The stochastic SRIF time 
update uses Givens transformations which exploit the upper 
triangular structure of the R matrix. 

Both the UD and SRIF time updates use subroutines of the 
ESL for most of the matrix operations. The SRIF stochastic 
update includes calculation of the smooth gain arrays which 
must be saved if smoothing will be done later. For the UD fil- 
ter, these smooth gains are calculated in a separate routine 
which can be skipped if smoothing will not be needed. 


The OASIS/GIPSY filter has a number of options for eval- 
uating filter models used for parameter estimation and covari- 
ance analysis. These include consider analysis tor both the UD 
and SRIF formulations, and a unique evaluation mode avail- 
able only with the UD filter. 

Consider parameters are usually bias parameters not included 
in the estimated parameter state vector. From the measure- 
ment partials, the sensitivity of estimated parameters to these 
considered (not estimated) parameters can be computed in the 
filter. There are several reasons for not estimating a parameter: 
(1) certain parameters, such as fiducial station locations, may 
be held fixed in order to define a reference frame and/or 
length scale; (2) it may be computationally impossible to adjust 
certain parameters, such as all the coefficients in a gravity 
field; (3) a physical effect cannot be modeled adequately to 
produce a reliable estimate, yet it is still desirable to calculate 
the penalty for leaving it out of the filter model (state vector). 
By including the effects of considered parameters, the ordi- 
narily overly optimistic computed covariance for the estimated 
parameters is degraded somewhat depending on the sensitivi- 
ties, and a more realistic covariance may result. 

Let y be the unestimated considered parameters with mea- 
surement partials A^. Then the sensitivity is [14] 

S = = _p A Tp-i A ( 20) 

dy c v c 
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where the ^ denotes the quantities computed from the filter 
measurement and time updates without including the effects 
of the y c parameters, and ? v is the measurement noise covari- 
ance matrix. The consider covariance , which includes effects 
from both estimated and considered parameters, is defined as 

P = P + SP S T (21) 

con c x ' 

where P c is the covariance matrix for the considered (unad- 
justed) parameters. 

Both the SRIF and UD filter implementations^ the GIPSY/ 
OASIS software can handle any parameter as estimated or con- 
sidered. However, only the UD filter has more general filter 
error evaluation capabilities, of which consider analysis is one 
type. The general error evaluation algorithm [24] permits 
evaluation of the effects of designing a filter without the cor- 
rect data noise, a priori covariance, data weights, or process 
noise model. A truth model and filter model are specified, and 
the evaluation mode of the UD filter allows the user to assess 
the effects of mismodeling on filter accuracy. Under real con- 
ditions, the truth model is not actually known, so a filter 
analyst usually tries to look at a range of possible truth models 
to see if, over a range of reasonable truth models, the filter 
results are especially sensitive to any parts of the filtering 
strategy being used. The goal is to design a filter model which 
is reasonably accurate based on available information about 
the problem, but which is also stable so that the results will 
not be strongly degraded if slightly incorrect filtering param- 
eters are used. In [25] is an example of a comprehensive filter 
evaluation which studied the use of stochastic solar radiation 
pressure models for GPS orbit determination. A number of 
process noise inputs were identified for which the filter 
showed little sensitivity to mismodeling. These process noise 
models were later used successfully for estimation of small 
GPS accelerations over arcs of several weeks (also see Sec- 
tion IV of this article). 

In the evaluation mode, the filter uses suboptimal Kalman 
gains saved in an evaluation file from an earlier filter which is 
run purposely with what is believed to be an incorrect model 
in order to generate suboptimal gains. The measurement up- 
date is expressed as 

P = (I - KA)P + a(K - K)(K - K) T (22) 

where the optimal Kalman gain is K and the arbitrary gain 
from the evaluation file is K. Let P = (I - KA)P and UDU T = 
P. Then K and a are computed using the measurement update 
formulas in Eqs. (2-4), and the vector X = K - K is formed. A 
rank-1 update computes the UD covariance factors: 


UdO t = UDU t + aXX T (23) 

The time update in the evaluation mode follows the same 
form as the original filter time update except that the original 
filter stochastic time constants and process noise sigmas are 
replaced with the evaluation mode time constants and process 
noise sigmas. 

F. Performance Comparisons Between the UD Filter 
and the SRIF 

The dual-filter option in the OASIS/GIPSY software facili- 
tates comparisons between the UD and SRIF programs. Many 
tests have been conducted to examine the numerical stability 
and speed of these two filters. At this time, they appear to be 
equally stable and numerically sound. However, there can be a 
substantial difference in processing speed, depending on the 
type of problem being filtered. Our experience with timing the 
UD and SRIF algorithms matches very closely the predictions 
in [14] . The SRIF measurement processor is about 33 percent 
faster than the UD measurement processor. For time updates, 
however, the UD algorithm is substantially faster when smooth- 
ing gains are not needed. If smoothing gains are computed in 
the filter for later use in smoothing, then once again the SRIF 
runs faster. Smoothing itself takes essentially the same amount 
of time whether the UD or SRIF filter was run first. For any 
specific filter run, the comparison between SRIF and UD run 
times depends on the number of stochastic parameters, the 
number of batch intervals, and the number of measurements. 
For typical GPS orbit estimation scenarios involving the geo- 
detic and low-earth orbiter applications described below, the 
SRIF is as much as 50 percent faster than the UD filter for 
complicated cases requiring smoothing. If the sensitivities in a 
consider analysis are needed for stochastically estimated 
parameters, the sensitivity matrix must also be smoothed. 
Sensitivity smoothing in the OASIS/GIPSY filter occurs in a 
natural and efficient way with the SRIF, while the UD sensi- 
tivity smoothing is usually substantially slower. The UD filter 
is about 25 percent faster than the SRIF in simple runs when 
no smoothing is needed. 

With the UD formulation in OASIS/GIPSY, the analyst is 
permitted to have the variance of a parameter equal zero. In 
certain analyses, such as a case where it is desired to constrain 
a parameter to a specific value, this capability can be useful. 
On the other hand, the SRIF formulation permits zeros on the 
diagonals of the information array, corresponding to an infi- 
nitely large covariance. As discussed above, the SRIF tends to 
be faster in situations when many measurements are being 
processed or when smoothing is needed. Only the UD filter, 
however, has an evaluation mode for studying filter mismodel- 
ing. The decision to fully implement both filter mechaniza- 
tions in the OASIS/GIPSY software allows the analyst full 
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flexibility in the selection of the most efficient and appro- 
priate filter for each specific situation. 

IV. GPS Positioning Applications 

Even though only seven GPS satellites are currently operat- 
ing, high-accuracy geodetic results have already been obtained. 
In addition to the geodetic studies carried out with the exist- 
ing limited set of satellites, numerous covariance analyses and 
simulations have been conducted in preparation for the full 
21 -satellite constellation expected to be operational in the 
early 1990s. A number of new approaches to high-accuracy 
earth-orbiter and ground positioning applications show great 
promise for demanding projects of the future such as TOPEX 
and Eos. In this section, a sampling of these results and new 
tracking techniques will be presented, with emphasis on how 
the orbit determination and filtering strategy play a role in 
GPS-based positioning and navigation. 

A. Ground Station Positioning Results 

A number of GPS campaigns and experiments took place 
in 1984 and 1985 [26, 27], in which typically 5-10 GPS 
receivers were located at sites of geodetic or geophysical inter- 
est. Baseline results showing precision and accuracy of 2-4 
parts in 10 7 relative to baseline length were initially reported 
[27-29] , with quality of the baselines generally assessed by 
comparison to VLBI and SLR, and by daily repeatability of 
the solutions. To reach these accuracies, GPS orbits were 
improved through estimation, generally with single-batch 
least-squares fits with doubly differenced dual-frequency 
carrier phase. Double differencing is a technique for cancella- 
tion of GPS and station clock errors which requires simulta- 
neous mutual visibility. In most, but not all of these cases, 
fiducial or reference ground sites were held fixed where GPS 
receivers were collocated with very precisely determined VLBI 
antennas in order to establish a reference frame and length 
scale. Without orbit improvement, GPS positioning results 
obtained using broadcast ephemerides only were, in most 
cases, about an order of magnitude worse. 

In [30] , results were presented from one of the 1985 GPS 
experiments showing accuracy and daily repeatability of 1-2 
parts in 10 7 for baselines up to 313 km. This improvement 
was due to a number of factors, including combining more 
than one eight-hour daily pass for orbit estimation, bias fixing, 
and an improved fiducial network geometry. Bias fixing is a 
technique for applying the constraint that the carrier phase 
range ambiguity be an integer number of wavelengths. This 
ambiguity is determined for as many baseline-satellite com- 
binations as possible, and the overall solutions are readjusted 
while the ambiguities are fixed at their integer values [31 , 32] . 
One approach to bias fixing [33] uses the SRIF formulation 


to optimally adjust all the biases and other estimated param- 
eters and to update the covariance matrix with an algorithm 
which is fast and does not require iterating or expensive re- 
filtering. The major weakness with solutions from only carrier 
phase data is a relatively low accuracy in the eastern baseline 
components due to the predominantly north-south GPS satel- 
lite tracks. Bias fixing can largely remove this weakness in 
the eastern directions. 

Kalman (factored covariance) filter-oriented approaches to 
the GPS estimation problem with emphasis on orbit determi- 
nation further improved baseline results to the level of 2-4 
parts in 10 8 [34] for baseline distances of more than 1000 km. 
The key strategies emphasized included multi-day (one-week) 
Kalman filtered orbit solutions, simultaneous parameter esti- 
mation, determination of three (instead of the usual two) GPS 
solar radiation pressure parameters, random-walk models for 
zenith troposphere delay fluctuations, and simultaneous 
processing of carrier phase and pseudorange data. The pseudo- 
range data tightly constrain the clocks and carrier phase ambi- 
guities, improving the results particularly in the east for base- 
lines and down-track for orbits. The transmitter and receiver 
clocks are estimated as white noise parameters, leading to 
essentially the same results as double differencing but with 
considerably less complexity. Direct comparisons between 
independent GPS orbit solutions [34] indicated orbit preci- 
sions of 1-2 m had been achieved. Further refinements to the 
filtering strategies have produced improved GPS accuracies of 
better than 1 m for orbits and 1-2 parts in 10 8 for baselines 
up to 2000 km [35]. These additional refinements included 
modeling of a GPS maneuver which enabled the data arcs to 
be extended from one to two weeks and constrained stochastic 
GPS force parameter estimation. Figure 5 shows locations of 
ground receivers from North American GPS experiments in 
1985 and 1986. Figure 6 shows the improvement in daily 
baseline repeatability using two-week orbit arcs when process 
noise models for force parameters were used. For arcs of one 
week or less, it was not necessary to use stochastic force 
models. In these experiments, all the measurements were con- 
fined to the same eight-hour interval each day and it was not 
possible to determine whether stochastic solar pressure or 
stochastic three-dimensional thrust parameters fit the data 
better. It is expected that in the future, with longer arcs and a 
more global tracking network, the physical nature of these 
stochastic forces will be better understood. 

Orbit verification to better than 1 m is illustrated in Figs. 7 
and 8. In each test, the data were partitioned into one -week 
arcs and independent orbit solutions were obtained and com- 
pared. The interleaved arc orbit comparison (Figs. 7a and 
7b) shows sub-meter repeatability for two well-tracked satel- 
lites. A more strenuous test using orbit prediction is shown 
in Figs. 8(a) and 8(b), demonstrating sub-meter agreement for 
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GPS 8. Figure 8(b) also shows how the results improve when 
pseudorange is processed with carrier phase and stochastic 
models for tropospheric delay fluctuations are used. 

A convincing test of GPS-based positioning accuracy comes 
from comparisons between 2000-km GPS-determined base- 
lines with VLBI measurements of the same baselines. VLB1 
is a completely independent geodetic system referenced to the 
inertial quasar -defined frame, whereas the GPS-based tech- 
niques use fiducial points tied to the VLBI system through 
local surveys to determine accurate GPS orbits, thereby trans- 
ferring fiducial control to other stations whose coordinates 
are estimated. Figure 9 shows that 2000-km baselines deter- 
mined with GPS agree with VLBI to better than 1 .5 parts in 
10 8 accuracy, which is close to the advertised accuracy of the 
VLBI system itself. 

Recent improvements in bias fixing techniques [33, 36, 
37] have enabled successful ambiguity resolution over regions 
up to 2000 km in extent. In these cases, cm-level agreement 
between GPS and VLBI has been achieved in horizontal base- 
line components with single-day arcs, provided that the sta- 
tions are spaced appropriately in these regions so that ambi- 
guity resolution can proceed over all or nearly all station- 
satellite combinations. To reach high accuracy in the vertical 
baseline components as well, multi-day arcs seem to be needed 
at the present time, but with worldwide tracking networks 
and additional GPS satellites, shorter tracking arcs will be 
adequate. GPS solutions seem to show lower scatter than 
VLBI solutions in the vertical component, perhaps because 
the GPS measurements are robust enough to estimate the 
troposphere stochastically with sub-centimeter precision. 

B. Earth Orbiter Positioning with GPS 

The TOPEX/Poseidon satellite [38] , scheduled for launch 
in late 1991, will lead to significant progress in the study of 
the interaction of the oceans and climate on a global scale. To 
understand global weather patterns, experiments like TOPEX 
are needed to study the world ocean’s circulation. The oceans 
redistribute heat, with warm currents carrying one-half of the 
excess heat from the tropics to the poles. The seas hold much 
more heat than the atmosphere and moderate the seasonal 
temperature fluctuations. Global circulation is observable from 
space because ocean movement causes bulges and depressions. 
The sea surface height can be measured with satellite altim- 
eters, as Seasat demonstrated in 1978. The detailed sea surface 
topography is a complicated function of ocean currents and 
the geopotential. A complete analysis of TOPEX data will 
separate the mean sea surface, or marine geoid, from time- 
varying currents, using repeat orbits and averaging techniques. 
The ability to map global ocean currents with TOPEX will 
have far-ranging benefits affecting weather and climate predic- 


tion, fishing, commerce, and shipping. The El Nino effect, 
which is a huge atmospheric seesaw leading to weather rever- 
sals, floods, economic hardship, heavy rain, and droughts, may 
be predictable in advance and eventually better understood 
through monitoring of the ocean surface with satellites like 
TOPEX, since it is typically accompanied by a raised mean 
ocean level of about 20 cm and an increase in the ocean tem- 
perature of about 2 deg C off the western coast of South 
America. 

TOPEX will fly in low-earth orbit at about a 1334-km alti- 
tude, carry an altimeter precise to about 2 cm to determine 
the range between the satellite and the sea surface, and will 
map the topography of the ocean surface. A key to the success 
of the mission is that the TOPEX orbit altitude error be kept 
below about 13 cm during the mission. Several different orbit 
determination techniques will be used on TOPEX, including 
ground-based laser ranging, but TOPEX will also carry a GPS 
flight receiver as part of an experiment to demonstrate precise 
low-earth satellite orbit determination. The TOPEX GPS demo 
is the first of the anticipated high-precision earth orbiter GPS 
applications for satellites at a wide range of altitudes, including 
geosynchronous as well as highly elliptical orbits. 

The major error source for differential GPS dynamic track- 
ing of TOPEX [1] is uncertainty in the model for the geo- 
potential. This is also expected to be the limiting error for 
other tracking techniques [38] . Recent attention has focused 
on an orbit determination approach using differential GPS in 
which the orbit filter is designed to desensitize the results as 
much as possible from gravity errors [39] . A three-dimensional 
fictitious force is modeled as process noise and estimated in 
this approach, referred to as reduced dynamic tracking. The 
role these stochastic force parameters play in the orbit filter is 
controlled through the process noise r and a w (Eqs. 11-14), 
and depends on the accuracy of the dynamic models for the 
satellite. If the dynamic models are very poor, the proper 
strategy is to set r 0 and o ss °°, and this is the limiting 
case of non-dynamic tracking. The other limiting case is 
purely dynamic tracking in which maximal weight is given to 
the dynamic models and the stochastic force parameters are 
essentially turned off (r -> 00 and o ss = 0). In between is 
reduced dynamic tracking , where a finite r and o ss determine 
how much reliance is placed on the dynamics and to what 
extent the unmodeled forces will be fitted out with the data. 
Although gravity is the major error source for TOPEX, the 
techniques used to minimize gravity-related orbit errors on 
TOPEX will also be useful for minimizing virtually any mis- 
modeled force for other satellites using differential GPS 
tracking. 

Figure 10 compares predicted performance of non-dynamic, 
reduced-dynamic, and dynamic techniques for TOPEX orbit 
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determination with GPS. Six globally distributed GPS ground 
receivers and the one TOPEX flight receiver were assumed to 
produce 5-cm pseudorange and 0.5 -cm carrier phase over five- 
minute data intervals in a single two-hour arc for this covari- 
ance study. This level of performance can be achieved only 
with advanced receiver and antenna design, but current results 
with the Rogue receiver/antenna are very encouraging and 
indicate that these accuracies can be met (Fig. 3). The pre- 
dicted altitude error shown in Fig. 10 for the dynamic solution 
strategy is dominated by the gravity field error, assumed to be 
50 percent of the difference between the GEML2 and GEM 10 
fields. A comparison with the covariance matrix for the 
GEM-T1 1 40 1 gravity field shows that representing the gravity 
error with the 50-percent geopotential difference for this 
analysis produces about the same results as produced with the 
current (GEM-Tl ) gravity field covariance. The large (~25-cm) 
excursions for the dynamic tracking approach occur over 
regions of the earth where the gravity field is poorly known, 
such as over the oceans. However, these oceanic regions are the 
key regions for an oceanographic experiment such as TOPEX. 
The non-dynamic approach also suffers from occasional sud- 
den increases in the altitude error, primarily due to times of 
weaker viewing geometry. Since the non-dynamic tracking 
solutions are more data-intensive, they are more sensitive to 
times when viewing geometry is compromised. The reduced 
dynamic strategy performs best overall. The stochastic force 
parameter estimation (r, a ) was chosen to be optimal for a 
gravity field accurate to about the level of the current GEM-Tl 
field, but as pointed out in [39] . even if the accuracy of the 
gravity field is only approximately known, the results are 
insensitive to the stochastic model as long as a relatively con- 
servative approach is taken. This is one of the main advantages 
of the reduced dynamic tracking strategy. 

The reduced and non-dynamic techniques compensate not 
only for gravity mismodeling. but also for any type of force 
or acceleration affecting spacecraft motion, such as drag, 
radiation pressure, gas leaks, maneuvers, etc. Estimation of 
fictitious thrust parameters for the GPS orbits themselves 
when tracking over relatively long (two-week) arcs (Fig. 6) is 
a simple example of the application of reduced dynamic con- 
cepts to satellites in high-earth orbits. It was found that the 
results in Fig. 6 showed very little sensitivity to r and a ss over 
a wide range of values, as predicted in [39] . The GPS orbit 
strategy, however, still places very high weight on dynamic 
models since the constrained stochastic forces being estimated 
represent a tiny perturbation to the purely dynamic orbit 
determination. 

C. Other Related GPS Applications 

Certain GPS applications show great promise but require a 
global tracking network and a full GPS constellation, such as 
earth orientation monitoring [4] . Ultimately, the filtering and 


estimation strategy may be a major factor in the utility of GPS 
for earth orientation because one of the advantages of GPS is 
the high time resolution it provides. Another example of the 
importance of Kalman filtering is the determination of tropo- 
spheric zenith delay parameters as part of the GPS orbit 
determination process. GPS-based techniques have the poten- 
tial of tracking tropospheric fluctuations at sub-centimeter 
levels with resolution of a few minutes. A GPS-based network 
for earth orientation monitoring, troposphere calibration, and 
determination of geocentric station coordinates is being 
studied for use in the Deep Space Network since these are 
major limiting errors for deep space navigation and planetary 
exploration. GPS tracking will be able to determine the loca- 
tion of the geocenter relative to any set of ground stations 
with an accuracy of better than 5 cm [5] , in addition to its 
cm-level relative station positioning capability. The NSWC 
precise ephemeris. which is determined from a worldwide 
Air Force tracking network and is more accurate than the 10- 
20-m ephemeris broadcast down on the P-code. is currently 
determined along with values for earth orientation parameters 
[41] using a SRIF and RTS smoother [42] and process noise 
models for polynomial clock coefficients, troposphere param- 
eters. and solar pressure coefficients. GPS also has the poten- 
tial for global time transfer at the sub-nanosec level, which 
would hold great benefits for the Deep Space Network and 
other radio astronomy observatories. 

If GPS antennas are placed at different points on an air- 
plane, the attitude of the aircraft can in principle be deter- 
mined from GPS measurements. An experiment was recently 
conducted [43] to demonstrate aircraft position, velocity, and 
attitude determination using GPS for a synthetic aperture 
radar experiment to study ocean currents, and the data are 
being analyzed at JPL. GPS receivers can also be placed on 
ships and buoys for research in seafloor geodesy [44] , a rela- 
tively new field. 

There are numerous real or near-real-time non-military GPS 
applications in which the cm-level accuracies required for 
geodesy or low-earth orbiter positioning are not essential but 
for which estimation speed is important. The NASA Sympo- 
sium on GPS Space Applications, held in Pasadena, California, 
November 18-19, 1987, emphasized many of the near-real- 
time positioning and navigation GPS applications for space- 
craft maneuvering near the Space Station and for other earth 
orbiters, including those at high earth altitude. Since many 
earth-orbiting satellites launched in the future will carry GPS 
equipment, including the Space Shuttle, it is likely that GPS 
will have an increasingly visible role in a variety of different 
missions. Although observing GPS from high earth altitude 
poses some special visibility problems, detailed analysis pre- 
dicts that meter-level orbit accuracy can be achieved even for 
geosynchronous satellites [6] . 
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In the future, it is expected that factorized filtering tech- 
niques will be applied to problems which would strain present 
computing resources to the limit. One example is recovery of 
gravity coefficients from GPS and low earth orbiter differen- 
tial tracking data. Because there can be thousands of gravity 
coefficients and the data arcs required include many satellite 
orbit revolutions, the number of estimated parameters is 
expected to be in the tens of thousands [45] . Numerical 
stability in these situations demands extremely stable, factor- 
ized filtering algorithms. Large problems such as this one 
would be intractable, even on a supercomputer, were it not 
possible to take advantage of the local block structure in the 
factored covariance matrix and the fact that much of the 
matrix is filled with zeros [45] . 


V. Summary 

The GPS multi-parameter estimation problem with orbit 
adjustment is well suited to a linearized batch sequential filter. 
Both the SRIF and UD-factorized filters have been used for 
high-accuracy GPS applications. Although neither is preferred 
over the other with regard to numerics, for different situations 
one or the other is sometimes more convenient or less expen- 
sive from a computational viewpoint. Determination of 
station locations to the cm level is currently feasible using 
GPS, with accuracy rivaling other state-of-the-art geodetic 
techniques. Accuracies of 1-2 parts in 10 8 relative to base- 
line lengths up to 2000 km have been demonstrated over base- 


lines measured independently by radio astronomy interfero- 
metric methods, and GPS orbits determined simultaneously 
can now be estimated to better than 1-m accuracy. Key as- 
pects of the estimation strategy include: collocation of three 
or four GPS receivers at fiducial sites with a priori, well-known 
coordinates; process noise modeling of fluctuating quantities 
such as clocks, tropospheric delays, and small but significant 
unmodeled satellite accelerations; multi-day orbit solutions; 
combined processing of pseudorange and carrier phase; and 
bias fixing techniques to apply integer constraints to carrier 
phase ambiguities. Future high-accuracy GPS applications 
also include sub-decimeter orbit determination for satellites 
such as TOPEX and Eos with non-dynamic and reduced 
dynamic tracking techniques, which use process noise to de- 
sensitize the filter to unmodeled or mismodeled forces. 

Although considerable progress has been made in just a few 
years in the development of strategies for precise positioning 
with GPS, improvements in GPS technology promise even 
more gains in the near future as advanced receivers and anten- 
nas become available. Global tracking networks will greatly 
increase the scientific potential of GPS, and will also lead to 
a dramatic increase in the number of estimated state param- 
eters. Precise pseudorange with noise of several cm will signifi- 
cantly advance virtually all the GPS positioning techniques 
developed so far. Such precise data and the anticipated filter- 
ing with very large, sparse matrices will probably require re- 
examination and further enhancement of the filtering ap- 
proaches currently being used. 
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L1{t) = P(t) B coslu^t) + C(t) B' sin(cj-|t) 


{LI CARRIER = 1.57542 GHz) 


Fig. 1 . GPS carrier phase is modulated by the P-code [P (t)] and the 
C/A code [C(t)]. The P-code contains nominal values for GPS orbits 
and clocks. Carrier phase and P-code pseudorange are used in 
high-accuracy applications. 


- (PSEUDORANDOM SQUAREWAVE CODE AT 1 023 Mbs) 




e USER POSITION MEASURED DIRECTLY 
WITH RESPECT TO GPS SATELLITES 

# MAJOR ERRORS 

GPS POSITION AND TIMING 
MEASUREMENT NOISE 

e INSTANTANEOUS POSITION ERROR 
C/ A CODE: 50-100 m 
P-CODE: 10-15 m 



e USER POSITION MEASURED WITH 
RESPECT TO PRECISELY KNOWN 
GROUND NETWORK 

# SUBSTANTIAL CANCELLATION OF 
GPS ORBIT AND CLOCK ERRORS 

• INSTANTANEOUS POSITION ERROR 

C/A CODE: 5-15 m 
P-CODE: 1 -3 m 


Fig. 2. (a) Direct GPS tracking, and (b) differential GPS tracking. 
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LINEAR POST-FIT RESIDUALS, 



0 1 2 


TIME, hr 

Fig. 3. Comparison between performance with different antenna and receiver designs: (a) post-fit residuals, mostly multipath, for dual- 
frequency pseudorange averaged over six minutes from a common T1 antenna/ receiver combination, and (b) pseudorange scatter obtained 
in January 1988 with the JPL Rogue receiver and a choke ring antenna designed to minimize multipath, for single-frequency data averaged 
over 2- and 30-minute intervals. 
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SRIF CPU TIME, mm 



BUFFER SIZE, number of measurements 

Fig. 4. Comparison of CPU time required for the SRIF for different 
measurement buffer sizes. 




EAST NORTH VERTICAL LENGTH 


Fig. 6. Daily repeatability for the 1314-km Mojave, California - Fort 
Davis, Texas baseline, for a two-week data arc in November, 1985. 
Lowest scatter is obtained with process noise force models. 


18 






13 14 15 16 17 18 19 



RMS DIFFERENCE BETWEEN 
ORBITS FROM TWO INTERLEAVED 
AND INDEPENDENT ARCS 



GPS 8 GPS 6 


13 14 15 16 17 18 19 20 21 22 23 24 



RMS DIFFERENCE 
BETWEEN ARC 1 
AND ARC-2 ORBITS 



CONSTANT CONSTANT STOCHASTIC 

TROP TROP TROP 


Fig. 7. (a) RMS orbit difference is computed from interleaved orbits 
determined independently from November 13-15-18 and Novem- 
ber 14-16-19; (b) RMS for GPS 6 and GPS 8 is shown computed 
over a 24-hour period on November 1 7 when no measurements were 
used for either solution as shown in (a). 


Fig. 8. (a) One week arcs in November 1985 used for orbit predic- 
tion and RMS computation. Orbits from the first week are mapped 
ahead and compared to orbits determined independently using data 
only from the second week; (b) prediction test (Fig. 8a) for GPS 8 
shows RMS well below 1 meter in all three components, with best 
results obtained using stochastic troposphere models and com- 
bined carrier phase and P-code pseudorange. RMS is computed 
over a 6-hour interval. 
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Fig. 10. Predicted altitude accuracy for TOPEX orbit determina- 
tion using dynamic, reduced dynamic, and non-dynamic tracking 
techniques. 



FIDUCIALS. HAYSTACK, FORT DAVIS, HAT CREEK 


Fig. 9. Comparison between GPS and VLBI independent mea- 
surements of 2000-km baselines shows agreement of 1 -2 parts 
In 10® in all vector components: (a) Hat Creek- Fort Davis, and 
(b) Richmond -Haystack. 
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